Transverse wavevector dependent and frequency dependent 
dielectric function, magnetic permittivity and generalized 
conductivity of interaction site fluids. MD calculations 

for the TIP4P water 

Igor P. Omelyan 

Institute for Condensed Matter Physics, the National Ukrainian Academy of Sciences, 
1 Svientsitsky St., UA-290011 Lviv, Ukraine * 

Abstract 

It is shown that the dielectric properties of interaction site models of po- 
lar fluids can be investigated in computer experiment using not only the 
charge fluctuations but also correlations corresponding to a current of mov- 
ing charges. This current can be associated with a generalized dynami- 
cal polarization or separated into electric and magnetic components. The 
first approach deals with the dielectric permittivity related to a generalized 
conductivity, whereas the second one leads to the functions describing po- 
larization and magnetization fluctuations separately. The last way is only 
the source to calculate the magnetic susceptibility for a system of interac- 
tion sites. The transverse wavevector- and frequency-dependent dielectric 
functions and magnetic susceptibility are evaluated for the TIP4P water 
model in a very wide scale of wavelengths and frequencies using molecular 
dynamics simulations. We demonstrate that the transverse part of the di- 
electric functions may differ drastically from their longitudinal component. 
A relationship between the two approaches is discussed and the limiting 
transition to the static dielectric constant in the infinite-wavelength regime 
is analyzed. The propagation of transverse electromagnetic waves in the 
TIP4P water is also considered. 
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I. MOTIVATION 



In a recent paper [1], a computer adapted fluctuation formula suitable for self- 
consistent calculations of the longitudinal wavevector- and frequency-dependent dielectric 
function e L (k, u) for interaction site models (ISMs) of polar systems has been proposed. 
As a result, a detailed analysis of this function in the entire wavelength and time scales 
was carried out for the TIP4P water model using molecular dynamics (MD). It was shown 
that the choice of microscopic variables for the operator P of polarization density plays an 
important role in a correct reproduction of dielectric properties, and only true microscopic 
variables, which explicitly take into account the charge distribution within molecules, can 
be used to determine the frequency dependence of the dielectric permittivity of ISMs at 
arbitrary wavenumbers. 

It is a common practice to investigate the dielectric properties of ISMs on the basis of 
charge fluctuations [1-3]. In the infinite- wavelength limit k — > 0, these fluctuations reduce 
to the well-known longitudinal dipole moment correlations, which are usually considered 
in computer experiment to obtain £ L (k, u) at zero and small wavevector values [4-12]. The 
transverse dipole moment fluctuations are used sometimes [5, 9, 10-12] for treating the 
transverse component e T (k,u) of the dielectric permittivity. However, such an approach, 
being exact for point dipole systems [13, 14], cannot be applied to calculations of the 
genuine transverse dielectric function of ISMs at nonzero wavenumbers. 

The problem of computations of e T (k, u) for ISM fluids is due to difficulties [15, 16] 
when constructing the transverse part of P. While the longitudinal part Pl can easily 
be expressed in terms of the operator of charge density, to define the transverse compo- 
nent Pt involving additional dynamical variables is necessary. There are two approaches 
for describing the electromagnetic phenomena in ISMs. They differ between themselves 
in the way of how to treat the current of moving charges. In the abbreviated approach 
[17], electric and magnetic parts of the total current are undistinguishable from one an- 
other in the presence of spatially inhomogeneous fields, i.e, when k ^ 0. In such a 
case, the magnetic part is included into a generalized P-vector and the electromagnetic 
phenomena in the system are determined by a transverse dielectric permittivity e T (k,u) 
which is directly connected with the generalized conductivity a T (k, uj) = ^(e T (k, u>) — 1). 
The second approach [18] is based on a separation of the microscopic current into elec- 
tric and magnetic parts at arbitrary wavenumbers. This leads to the introduction of 
two functions, e T (k,u) ^ e T (k,u) and fi(k,u), describing the transverse fluctuations 
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of polarization and magnetization densities separately. The last function is associated 
with a generalized magnetic permittivity of ISM systems. In this approach the trans- 
verse dielectric permittivity e T (k,u>) can easily be reproduced using functions e T (k,u) 
and n(k,u). In the infinite- wavelength limit k — > 0, when the spatial dispersion can 
be neglected, the two approaches become completely equivalent. In this case the elec- 
tric phenomena are uniquely described by the frequency- dependent dielectric constant 
e{uj) = lim fe ^ e T (k, u) = lim fe ^ e L T (k, u), whereas a magnetic state is determined by the 
magnetic permittivity fi(tu) = lim fc _^ fJ>(k, u ) = [1 ~~ ^ nm fc^o(^ T (^) ^) ~ £ - L (k,uj))/k 2 ]~ 1 , 
where c denotes the velocity of light [17]. 

Until recently, quite a few papers [19, 20] dealt with the investigation of correlations 
related to the current of charges in ISMs. In these articles, spectra of the longitudinal 
and transverse components of the hydrogen current, on the time scale peculiar to the 
librational dynamics of water, were calculated for the TIP4P potential at small wavenum- 
bers. The first calculation of the transverse dielectric function e T (k) = lim^^o e T (k, u) 
has been performed by Raineri and Friedman [16], but in the static regime only and for 
the simplest £DS model. At the same time, there were no attempts to investigate the 
entire wavevector and frequency dependence of the transverse dielectric function and the 
magnetic susceptibility for ISM fluids. 

In the present paper the dielectric properties of ISMs are investigated on the basis 
of current correlations, including their separation into electric and magnetic parts. This 
allows us to determine both the longitudinal and transverse components of the dielectric 
permittivity e(k,uj) as well as the dielectric e T (k,u) and magnetic fi(k,uj) functions. 
Actual MD simulations are performed for the TIP4P model of water using the interaction 
site reaction field (ISRF) geometry [3] and the Ewald method. The results obtained for 
e T (k,u) and /i(k,ui) are presented in a wide region of wavevectors and frequencies. 

II. ELECTROMAGNETIC FLUCTUATION FORMULAS FOR ISM FLUIDS 

A. Generalized dielectric constant and conductivity 

We shall consider a polar fluid consisting of N identical molecules which are composed 
of M interaction sites and enclosed in a volume V. The microscopic density of charges for 
such a system at point r G V and time t is of the form Q(r, t) = Y^iLi J2^Li Qa^( r — r i'(^))) 
where q a and r^(t) are the charge and position of site a within molecule i, respectively. 



3 



Shifting from the real coordinate space {r} into the {^-representation by the spatial 
Fourier transform {k} = f {r}e" ifc ' r dr, we obtain that Q(k,t) = J2?=i Ea=i q a e~ ik<(t) - 
For the investigation of dielectric properties, it is convenient to introduce the microscopic 
vector P(k,t) of polarization density. The longitudinal part P L = kP-k of this vector 
can be determined using the relation Q(r,t) = — divP(r,t) in fc-space, i.e., Q(k,t) = 
-ik-P(k,t). Then we find 

p L (fe, t) = %Q{k, «) = |ee ^~ ik - rm , (i) 

i=l a=l 

where k = k/k is the unit vector directed along k. 

Recently [1], it has been shown that the longitudinal component e L (k,u) of the 
wavevector- and frequency-dependent dielectric tensor can be calculated in computer ex- 
periment via the fluctuation formula 



e L (k,u)-l _ /L2L(-G L (M)) 



= ^(-<7 L (M)) . (2) 



e L (k,uj) l + hD(k)J?U-Gdk,t)) 

Here, G L (k,t) = ^P L (fe, 0)-P L (— k, t)^ jNd 2 is the longitudinal wavevector-dependent 

dynamical Kirkwood factor computed in simulation for a finite sample, ( ) Q denotes 

the statistical averaging in the absence of external fields, d designates the permanent 

magnitude of molecule's dipole moment dj = 1a r ii the Laplace transform is defined as 

~^L({0) = / {0 6"^d£, k-Q and T are the Boltzmann's constant and the temperature 
Jo 

of the system, respectively, h = AuNd 2 jVk^T and G^(k,t) = dGi J (k,t)/dt. The function 
D{k) takes into account boundary conditions, applied in simulation, and for finite samples 
it always is equal to 1 at k — 0. For nonzero wavevectors this function tends to zero in the 
case of Ewald summation, while D(k) = 3j 1 (kR)/(kR) within the reaction field geometry 
[1, 3], where R and j ± (z) = sm(z)/z 2 — cos(z)/z are the cut-off radius and spherical Bessel 
function of first order, respectively. For infinite systems (R — > oo) the function D(k) is 
equal to zero at arbitrary wavenumbers and the computer adapted formula (2) reduces 
to the well-known fluctuation formula in terms of the infinite-system Kirkwood factor 
g L (k,t) = lim A r_ 00 G' L (A;,t). 

As we can see, using the microscopic operator Q(k,t) of charge density allows one to 
determine uniquely only the longitudinal part of polarization vector P(k, t). To construct 
the transverse part, it is necessary to involve additional dynamical variables. We shall 
show now how to derive fluctuation formulas which give the possibility to computer both 
the longitudinal as well as the transverse component of the dielectric constant of ISMs. 
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The processes of dynamical polarization in dielectrics cause a current of charges. This 
current is described by the microscopic operator of current density I and may be associ- 
ated with the generalized polarization current [17], so that 

N,M i 

I(k, t) = £ ga V?(0e-*<W = -P(k, t) , (3) 

i, a 

where V"(i) denotes the velocity of site a within the molecule % at time t. The relation (3) 
can be considered as a more general definition of P-vector, because it allows to determine 

a --v yv y> a a a a a a 

its longitudinal P L = kP-k and transverse P T = (1 — kk)P = kx [Pxk] components, 
where 1 is the unit tensor of the second rank. It is worth to emphasize that the generalized 
vector P of polarization induction includes both electric as well as magnetic contributions 
(see a more detailed discussion of this point in the next subsection). In a particular case 
of longitudinal polarization this vector satisfies completely our previous definition (1) of 
Pl- Indeed, taking the full derivative of (1) with respect to time, it can be verified easily 
on the basis of equation (3) that ^P L (k,t) = kk-I(k,t) = I L (k,t). 

Let us apply to the system under consideration an external electric field E (k, t) which 
contains longitudinal Eq and transverse Eq components. The total field in the system 
E(k,t) consists of the external field and an internal field of charged sites. Neglecting 
the relativistic terms due to the dynamical magnetic field of moving charges, the internal 
field can be presented [1, 3] in the purely longitudinal form — Att(1 — D(k))P L . Therefore, 
the longitudinal and transverse components of the total field E = E L + E T are equal 
to E L (k,t) = Ev(k,t) - 4tt(1 - D(k))P L (k,t) and E T (k,t) = E^(k,t), respectively. 
The longitudinal e L (k,u) and transverse e T (k,u) components of the wavevector- and 
frequency-dependent dielectric tensor e{k,uj) = e T (k,u)l + (e L (k,u) — e T (k,u))kk can 
be defined via the material relations P LtT (k,u) = j^(e LT (k,u) — 1) E^T(k,u), where 
P LjT {k,uj) = ^P Li T(fe,<^)) and E L ^(k,uj) = ^-El,t(&, are macroscopic values of the 
polarization and total field, ( ) denotes the statistical averaging in the presence of the 

/oo 
{t}e~ lwt dt has been used for 
-oo 

functions P^^(k,t) and E^^(k,t). 

In the case of dynamical polarization, when I(k,to) = ^i(k,io)^ ^ at to ^ (note 
that static macroscopic currents of coupled charges are equal to zero), the Maxwell equa- 
tions can be written in a form as for conductors, namely, in terms of the generalized 
polarization conductivity a(k, uj). The longitudinal and transverse components <r L T (k, uj) 
of this wavevector- and frequency- dependent conductivity are defined via the material 
relations I^{k,uj) = a hl ,{k,uj)E^{k,uj). In the frequency representation, equation 
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(3) stays that i L ^(k,u) = icuP LtT (k, u). Thus we can express the dielectric constant in 
terms of the polarization conductivity as follows 

47T(7 T ryik.UO) 

e LT (*,a;)-l = . (4) 

' iuj 

According to the perturbation theory of first order with respect to external fields, we 
obtain for the macroscopic current J LiT (fc,cj) = Jzf icJ (^I L ^(k, 0)-Jl,t(— k, t)^ ^V^r ' 
where the multiplier {2} is to be included for the transverse component only. Then 
eliminating the electric fields and using relation (4) yields 

h<M^l - hsrucUKt)) h_ 

e h {k,uo) iuo + hD(k)3?UCL(k,t)) j^^W'^' ^ 
e T (k,u)-l = ^S? iu (.c T (k,t)) , (6) 

where 

(i L , T (fc,o).i L , T (-fc,o) 

c ^=- — — (7) 

are the longitudinal and transverse components of the wavevector-dependent dynamical 
Kirkwood factor of second order for the finite sample, whereas c hT (k,t) are the corre- 
sponding functions of the infinite system and c T (k,t) = C^ik^t). 

The fluctuation formulas (5) and (6) can be used in simulations to evaluate the lon- 
gitudinal and transverse components of the dielectric constant on the basis of equilib- 
rium current-current correlations (7). In the case of longitudinal fluctuations, the for- 
mula (5) is mathematically equivalent to the usual relation (2). Indeed, in view of the 
equality Ij,(k : t) = ^P L (fc,t), it can be shown easily that Ci,(k,t) = —^Gi,(k,t) = 
—Gi,(k,t). Then applying the Laplace transform we obtain J2L(Cl(A;, t)) = Gx(A;, 0) + 
Gx(fc, 0)) where Gi,(k, 0) = because the Kirkwood factor Gt,(k : t) is an even 
function on time, and we immediately recover the fluctuation formula (2) from (5). 
Furthermore, taking into account that Jzf iw (— G L (k, t)) = G L (k,0) — iuj^ iL0 (G L (k,t)) = 
•&iw(Ci,(k,t))/iu we obtain, in particular, that the static (t = 0) longitudinal Kirkwood 
factor G^(k) = G^k, 0) is connected with the dynamical Kirkwood factor of second order 
as 

G L (jfe) = l im ^UCdk,t)) = _ r° tC , L(M)dt ^ (g) 

IUJ Jo 

where the equality (8) holds for the infinite-system functions g L (k) and c L (k,t) as well. 
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Despite the difference, in general, between the longitudinal Kirkwood factors for the 
finite and infinite systems, we, nevertheless, can reproduce the infinite-system behaviour 
indirectly using the relations 

-2Lte L (M)) = J2L(G L (M)) ~ hD{k) ' ^(Mj) = J2L(C L (fc,t)) + hD{h) ' 

(9) 

Taking into account the Laplace boundary theorem lim w _ KX} iu;JSfi a ,(0(t)) = lim t _ +o 0(t), 
it can be shown from (9) that -^Gi,(k,i) ^ = -^g L (k,t) ^ and Cl(/c) = c L (k), re- 



t=o 



spectively. So that, contrary to the usual Kirkwood factor G^(k), the static Kirkwood 
factor of second order Cl(/c) is free of boundary effects. The transverse function c T (k,t) 
is equal to Ct(/c, t) at arbitrary times, because Et = Eq for nonrelativistic systems, and 
it can be obtained in simulations directly without additional manipulations. Moreover, 
the static Kirkwood factors of second order c LT (k) are presented in analytical forms (see 
Appendix). 



B. Transverse dielectric function and magnetic permittivity 

We now consider an alternative approach to describe the electromagnetic properties 
of ISM systems. This approach is based on the separation of macroscopic currents into 
electric and magnetic parts in the limit k,oo — > 0. Neglecting terms of order k 2 , uj 2 , ku 
and higher, the averaged values for the operator of current density (3) can be evaluated 
in (k, u;)-space as 

(I(k,u)) = ickx(M(w)) + iuo(V{uo)) , (10) 

where 'P(w) and jCi(u) are frequency components of the total electric dipole moment 
V(t) = J2iLidi(t) and the rotational part At(t) = J2iLi m i(t) of the magnetic dipole 
moment of the system, where rrti = j- Y^=\ Qa^i^^- I n our notations <5" = r" — and 
= V" — Vi = i7j x <5" are the positions and velocities of sites relatively to the molecular 
centre of mass and its velocity Vj, respectively, and 17, is the angular velocity of the 
ith molecule. 

Despite the fact that the separation (10) is realized uniquely only in the macroscopic 
regime at small wavenumbers and frequencies [17], we shall apply it at the microscopic 
level of description and extend to arbitrary k and uj. Such an extension can be performed 
by writing (in (fe, ^-representation) 
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I(M) =ickxM(k,t) + ^-V(k,t) , (11) 

where P(k,t) and A4(fc, t) are appropriate dynamical variables associated with the mi- 
croscopic polarization and magnetization densities, respectively. In view of equations (3), 
(10) and (11), these variables should satisfy the limiting transitions lim^ V(k, t) = V(t) 
and lim^o M.(k,t) = jVt(t). Moreover, taking the scalar product of equations (3) 
and (11) on unit vector k yields the condition T , i J (k,t) = P^(k,t), where the equal- 
ity k-[kx M] = has been used. 

The explicit expression for P(k,t) can be derived using the following procedure [16]. 
Let us split the operator Q(k,t) = J2iLi q i (k,t)e~ lk ' ri ^ of charge density into its molec- 
ular components q^k^t) = X^li q a Q~ lk ' S ^^- Then it is natural to introduce the total mi- 
croscopic polarization density as V(k,t) = J2iLiPi{k,t)e~ lk " ri( - t \ where p^k^t) is the po- 
larization density of the ith molecule, and to extend the relation Q(k, t) = — lk-P(k, t) to 
the molecular level, i.e., q.(k,t) = —ik'p^kjt). Further, taking into account the identity 
= 1 + £ Jq 1 due u ^, applied to the quantity £ = — ik-S®, and using the molecular charge 
electroneutrality E^li q a = °, we present ^(M) as ~ ik ' E^ii 9o*i(*) Jo due _iufc ' tf ? ( * ) . 
The last expression leads to Pi(k,t) = J2%Li ? o <^(0 Jo due~ iuk ' s i^ and, therefore, 

N M p -ik-8?(t) _ i 

?(M)-EE^) jr .a m e-^-^W . (12) 

i=l o=l 

Of course, such a procedure does not define P(k,t) uniquely, because q. is indifferent 
to the transverse part of p { . Nevertheless, it was assumed [16] to adopt equation (12) 
as the definition of microscopic polarization density for ISM fluids. From this definition 
it immediately follows that the longitudinal components of vectors P and "P coincide 
completely between themselves at arbitrary wavevectors. 

Since the functions I(k,t) and V(k,t) are already defined, equation (11) allows one 
to determine the transverse part of microscopic magnetization density, 

M T (M) = kx[M(k,t)xk] = J-(j(fc,i) - p(k,t))xk . (13) 

1C/C (XT 

It can be verified using equations (3), (12) and (13) that in the infinite-wavelength 
limit k — > the functions P(k,t) and Ai T (k,t) tend to the true microscopic variables 
J2? =1 di(t) = J2?=i Ea=i Q a Si(t) and kx 'E^ =1 [m i (t)xie] which correspond to the electric 
and rotational magnetic dipole moments of the system, respectively. However, at k ^ 
the vectors V(k,t) (as well as P(k,t)) and jCi T (k,t) take into account explicitly the 
charge distribution within a finite spatial extend of the molecule and, therefore, they can 
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not longer be associated with the point dipole densities &{k,t) = J2iLi d i (t)Q~ lh " ri ^ and 
^#r(fc,t) = kx Y$Li[™>i{t) xfc]e _lfc ' ri ^. For purely dipole models, when \q a \ — > oo and 
max a |(5"| — > 0, provided d — > const, we obtain P(k,t) — ► &(k,t), but J^(k,t) — > 
because of rrij — > and, thus, there is no magnetic response in this case, i.e., P{k,t) = 
P(k,t) at arbitrary wavevectors. 

In the present approach the electric phenomena are described by the frequency- 
dependent dielectric functions e LT (A;, u;), while the magnetic state is determined by the 
magnetic permittivity /i(k, uj) = /i T (/c, a;) (for the system under consideration the longitu- 
dinal magnetic susceptibility is absent). These quantities are defined via the material rela- 
tions P l ,tO, w) = j^(\ T (k, u)-l)E LtT (k, uj) and M T (k, uj) = ^(fi(k, uj)-l)H T (k, uj), 
where the statistical averaging of P = (P) and At = (Ad) is performed in the presence of 
external electric Eq and magnetic Hq fields. We mention that for nonrelativistic systems 
the total magnetic field H is indistinguishable from the external field, i.e., = Hq. 
We stress also that it is necessary to distinguish the generalized dielectric tensor e L T (k, uj) 
corresponding to the total current i" = dP/dt from the dielectric functions e LT (k, uj) re- 
lated to the electric part dV/dt of i". It is obvious that e L (k,uj) = e L (k, uj) because of 
V L (k,t) = P L (k,t), but e T (k,uj) ^ e T (k,uj) since, in general, P T (k,t) ^ P T (k,t). The 
fluctuation formula for the transverse dielectric function e T (k,uj) can be obtained in a 
similar way as for permittivities e LT (k,uj) (see the preceding subsection and Ref. [1]). 
The result is 

e T (fc, uj)-1 = h^(-g T (k, t)) , (14) 

where the infinite-system correlation function g T (k,t) = (p^(k, 0)-"Pt(— k, t)^ jlNd 2 
can be evaluated in simulation directly, because of E T = Eq. 

Finally, we would like to discuss about a relationship of the two approaches presented 
in more detail. In the abbreviated approach the electric and magnetic parts in the total 
current are not distinguished from one another. In such a situation the magnetic part is 
included into the generalized P-vector (3) and the transverse electromagnetic fluctuations 
are described by one function only, namely, by the generalized dielectric constant s T (k, uj). 
This including indeed can be realized because at k, uj ^ magnetic fields are expressed 
in terms of electric fields E using the Maxwell equation crot.E = —dB/dt in (k,uj)- 
representation, i.e., ckxE(k,uj) = —ujB T (k,uj), where B T (k,uj) = fi(k,uj)H T (k,uj). At 
the same time, the extended approach involves two quantities, e T (k,uj) and /i(k,uj), for 
the description. As far as these two approaches deal with the same microscopic current, 
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the quantities e T (k,u), e T (k,cu) and fj,(k,u) are not independent. The corresponding 
relation can be found comparing the right-hand sides of equations (3) and (11) between 
themselves and using material equations. Then one obtains 

^H^,^ — HL^- (15) 

and, therefore, the wavevector- and frequency-dependent magnetic permittivity is caused 
by the difference between functions e T (k,u) and e T (k,cu). 

The simplest way to obtain the explicit fluctuation formula for the magnetic permit- 
tivity is based on relation (15) and fluctuation formulas (6) and (14) for e T (k,u) and 
e T (k,u). In view of equation (11), the transverse current autocorrelation function c T (k,t) 
(7), appearing in fluctuation formula (6) for e T (k,u>), can be expressed in terms of the 
polarization g T (k,t) and magnetization s T (k,t) = (^A4t(&, 0)-A4t(— fc, £)^ q /2iV<i 2 cor- 
relation functions as c T (k,t) = c 2 k 2 s T (k,t) — d 2 g T (k,t)/dt 2 . It is necessary to point out 
that for spatially homogeneous systems, as in our case, the cross function (T't-M.t) , 
corresponding to correlations between the polarization and magnetization vectors, does 
not appear and it is equal to zero for arbitrary wavenumbers and times. Indeed, let us 
direct fe-vector along z-axis of the laboratory reference frame. Then it follows from the 
structure of equation (11) that Pt-AHt = \{V x M. y + V y M x ), where V x , y and M X)V are 
the xth. and yih components of vectors "P and jCi, respectively. Since the fluctuations 
of vector quantities in different directions of the fixed laboratory frame are statistically 
independent at equilibrium, we have that (^pT(k, 0)«A1t( — k, t)^ = 0. Further, as far 
as the relativistic effects have been neglected, the functions e T (k,u) and e T (k,u) are 
evaluated via fluctuation formulas (6) and (14), in fact, with a precision of order (f/c) 2 , 
where v c is the mean heat velocity of atoms. This leads to uncertainties of order 
c~ 4 in the evaluation of fj,(k,cu) via relation (15). Within the same precision we can putt 
(fj,(k,Lo) — l)/fj,(k,u) ~ fi(k,u) — 1, so that the desired fluctuation formula is 

n{k, uo)-l = -\uh^{s T {k, t)) + 0(cT 4 ) , (16) 

where the infinite-system function s T (k,t) can be reproduced directly in simulation, be- 
cause of H T = Hq. 

From the afore said, it is obvious that the two approaches are completely equivalent for 
the evaluations of the generalized dielectric permittivity e T (k, u) and they can be applied 
with equal successes to theoretical applications. However, the extended description is 
only the way to determine the magnetic permittivity at k ^ 0. It is interesting to 
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point out also that in the infinite-wavelength regime k — > 0, the functions e L (k,u) = 
e L (k,u>) and e T (k,u) differ between themselves by terms of order k A and higher, i.e., 
lim / t^o[^ L (^, ^0 - e T (^)^)]/^ 2 = 0. This statement can be examined on the basis of 
fluctuation formulas (2), (14) and explicit expression (12) for microscopic variable "P. 
Then, using relation (15), we obtain that in the abbreviated description the frequency- 
dependent magnetic permittivity fJ>(u) = lim^o fJ>(k, u) appears as a result of the limiting 
transition 

1 _ _J_ = ^ lim g T ( fc ^)- £ L ( fc ^) , 17) 
fi(u) C 2 k^O k 2 

which shows that differences between the longitudinal and transverse components of 
the generalized dielectric tensor e(k,u) in the infinite-wavelength limit are caused by 
magnetic properties of the system [17]. Thus, in the presence of spatially inhomo- 
geneous fields, the electromagnetic state is determined by the generalized transverse 
e T (k,u) and longitudinal e L (k,u) dielectric functions. When the spatial dispersion 
can be neglected, all electric and magnetic phenomena in the system are uniquely de- 
scribed by two quantities again, namely, by the frequency-dependent dielectric constant 
e{uj) = limfe^o £ L (k, u>) = lim^o e L T (k, uj) and magnetic permittivity [J,(u). 

III. NUMERICAL CALCULATIONS FOR THE TIP4P MODEL 

Molecular dynamics simulations were carried out for the TIP4P potential [21] at a 
density of 1 g/cm 3 and at a temperature of T = 293 K. We performed two runs cor- 
responding to the ISRF [3] and Ewald [22] geometries. The equations of motion were 
integrated on the basis of a matrix method [23] with a time step of 2 fs. The observation 
times over the equilibrium state were 500 000 and 1 000 000 time steps for the Ewald and 
ISRF geometry, respectively. The time correlation functions were calculated in an interval 
of 2 ps for the wavenumbers k — [0, 1, ... , 300]/c m i n , where fc m i n = 2ix / L = 0.319A 1 and 
L is the length of the simulation box edge. The parameters rj = 5.76/ L and /c max = 5/c m i n 
have been used in the Ewald summation of the Coulomb forces. Other details of the 
simulations are similar to those reported earlier [1]. 

A. Dielectric properties 

The longitudinal component G L (k) of the static wavevector-dependent Kirkwood fac- 
tor, calculated in the simulations within the Ewald and ISRF geometries, is shown 
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in fig. la by the full squares and dashed curve, respectively. Since, in the ISRF ge- 
ometry the function D(k) = 3j 1 (kR)/(kR) differs from zero considerably, especially 
at small wavenumbers, to evaluate the infinite-system Kirkwood factor g h (k) the per- 
formance of self-consistent transformations (9) is necessary. This result is plotted 
by the solid curve. Within the Ewald geometry [22] the function D(k) is equal to 
1 — /(f kj 1 (kp)(eric(r]p) + ^=pexp(—r] 2 p 2 )^Jdp — exp(— k 2 /Ar] 2 ), where the last term is 
to be included only if < k < k max , and at the given parameters of the summation D(k) 
is very close to zero (max fc ^ |-D(^)| < 0.00005). So that the infinite-system Kirkwood 
factor is equivalent to that, obtained in the simulations (excepting the case k = 0). As 
we can see from the figure, the both Ewald and ISRF methods lead to identical results 
for g L (k). The static Kirkwood factor of second order, Cl,t(^) (equation (7)), is pre- 
sented in fig. lb. As was pointed out earlier, this function is free of boundary conditions 
and the infinite-system dependence c LT (/c) can be reproduced directly in the simulations. 
This is confirmed by our calculations performed in the different geometries. The obtained 
values for Cl,t(^) are practically indistinguishable from those evaluated from analytical 
expressions (A5) for c LT (/c). They differ from one another within statistical noise only. 

Samples of the normalized, dynamical Kirkwood factor of second order, ty L:T (k,t) = 
Cl,t(&, 0/Cl,t(&), are plotted in fig. 2. The longitudinal infinite-system functions 
c L (k,t)/c L (k) within the ISRF geometry at arbitrary wavenumbers and for the Ewald 
geometry at k — have been determined applying the inverse Laplace transform to 
relations (9), whereas c T (k,t) = Ct(&, t) for the transverse functions. We note that 
C^(k,t) = Ct(/c, t) at k — 0. At k ^ the Ewald method reproduces directly the 
infinite-system behaviour. The agreement between the two sets of data for c LT (k,t), cor- 
responding to the ISRF and Ewald geometries, is quite good. It is worth to remark also 
that the zeroth time moment c L (k, t)dt is equal to zero for arbitrary wavenumbers 
(this statement directly follows from equality (8)), whereas, in general, J °° c T (k, t)dt ^ 
in the case of transverse total current fluctuations. 

As was shown in the preceding section, the longitudinal static Kirkwood factor G\,{k) 
of order zero can be determined through the first time moment on the dynamical Kirkwood 
factor C^(k,t) of second order (see equation (8)). The calculated in such a way function 
Gi,(k) within the ISRF geometry and the function g L (k) (via c L (k,t)) using the Ewald 
method are presented in fig. la by rotated and direct crosses, respectively. As we see from 
the calculations, these functions differ considerably from those obtained in the usual way. 
This difference occurs because the calculation of the expression Jzf iuJ (C L (k, t))/'\uj at small 
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frequencies is very sensitive to the precision of the evaluation of J^L(Cl(/c, t)) (dividing two 
small quantities between themselves). The problem is complicated additionally because 
this evaluation requires a numerical integration of time correlation functions which are de- 
fined in simulations approximately within a statistical noise. Therefore, this method is not 
recommended for the investigation of the longitudinal dielectric constant at low frequency 
values. However, in the opposite infinite-frequency limit, the calculation via equation (5) 
can be more convenient than using the usual way (2). For example, in this limit, where 
e h (k,u) — > 1, we obtain from (2) that J^ iui (—g L (k : t)) = g L (k) — iuJf iuJ (g L (k, t)) — > 0. The 
exact computation of iuJ% w (g L (k,t)) at large frequencies may lead to a problem (mul- 
tiplying of quantities with significantly different orders). This situation is absent in the 
case when fluctuation formula (5) is applied to calculations. 

And now we are in a position to discuss a behaviour of the generalized dielectric 
permittivity in the low frequency regime uj — > 0. As was established previously [1], the 
real part of the longitudinal dielectric permittivity e L (k,u) at uj — * tends to its static 
value £ L (k), while the imaginary part vanishes at arbitrary wavenumbers except the cases 
of two singularities, where e L (k) = ±oo. In the singularity range, the real part of the 
dielectric permittivity is equal to zero, whereas the imaginary part behaves as — i/ujrl° r (k), 
where r™ T (k) = J °° dt g L (k, t)/g L (k) is the longitudinal correlation time. Such a behaviour 
of the longitudinal dielectric constant causes the coefficient cr L (k) = lim^o <r L (k, uj) = 
l/A-KT^ik) 7^ for longitudinal conductivity in the singularity regions, whereas outside 
the singularities & L (k) = 0. The existence of a nonvanishing coefficient in the static limit 
does not lead, however, to macroscopic currents (as it must be for dielectrics), because 
the total longitudinal field in the system vanishes when |e L (/c)| — > oo. 

In the case of transverse fluctuations the pattern is qualitatively different. Accord- 
ing to fluctuation formula (6), the transverse dielectric permittivity e T (k,u) at low 
frequencies behaves as 1 + lim^o /i«5f iw (c T (A;, t))/'\uj = e'(k) — ie"(k), where e'(k) = 
3ft lim^o e T (k, uj) = 1 — h J °° tc T (k,t)dt denotes the real part, whereas the imaginary 
part e"(k) = — 3 lim^o £ T (k, uj) = Ana^ik) /uj is described by the wavevector-dependent 
generalized transverse conductivity a T (k) = lim^o -£iu(c T (k, t)) = c T (k, t)dt. 

We mention that the precision of calculations of the generalized dielectric permittivity 
at low frequencies is very sensitive to statistical uncertainties of data for total current 
fluctuations c LT (k,t) obtained in computer experiment. An additional source of errors is 
the truncation of long time tails in correlation functions at numerical integration. As a 
consequence, we could not provide a satisfactory reproduction of e T (k, uj) directly in terms 
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of c T (k, t), especially at small wavevectors and in the case of the Ewald method, where the 
length of the simulation was twice smaller than for the ISRF geometry. For the last rea- 
son, all other results will be presented in the ISRF geometry only. The computations show 
that a much more reliable evaluation can be performed when the generalized permittivity 
e T (k,uj) is found via relation (15), i.e., in terms of the transverse dielectric e T (k,u) and 
magnetic fi(k,u) functions which correspond to polarization g T (k,t) and magnetization 
s T (k,t) fluctuations. Then in view of (16), relation (15) transforms at u — > to 

2 1.2 

lim eJk, uj) = eJk) - ih- — lim Sf^sJk, t)) + 0(uj, c" 4 ) , (18) 

ui^Q UJ llI^O 1 

where e T (k) = lim^o e T (k, uj) = 1 + hg T (k) is the static dielectric function and 
g T (k) = lim t ^, 9 T (k, t). From equality (18) we immediately obtain that e' (k) = 
e T (k) - hc 2 k 2 / °° ts T (k, t)dt and e'^k) = Aita T {k)/uj with a T (k) = ±c 2 k 2 / °° s T (k, t)dt. 

The wavevect or- dependent dielectric function e T (k) is shown in fig. 3a by the solid 
curve. In the infinite-wavelength limit k — > 0, this function tends to the value e = 
lim k ^ e L (k) 53, corresponding to the usual dielectric constant in the absence of any 
spatial and time dispersions. It is equal to unity in the opposite limit k — > oo, indicating 
that the system exhibits no dielectric response with respect to strong spatially inhomo- 
geneous electric fields. The dielectric function e^ D (/c) = 1 + hg^ D (k), calculated in the 
point dipole (PD) approximation g^ D (k) = ( J 2 T (fc, 0)- J 2 ' T (— fc, 0)) /2A^<i 2 , is presented 
in fig. 3a by the dashed curve. As can be seen from the figure, this function behaves 
like that for a Stockmayer fluid [13]. The PD approach is suitable for the reproduc- 
tion of e T (k) at very small wavenumbers only, namely, at k <^ ir/r ~ 3.4A \ where 
r = max a |<5"| ~ 0.92A denotes the radius of the TIP4P molecule. In this wavevec- 
tor range V 2? and the spatial extend of charges within the molecule can be not 
taken into account when constructing the operator of microscopic polarization density. 
At greater wavevectors, the PD function differs drastically from that obtained within the 
exact interaction site description. In particular, the dielectric function e^ D (k) tends to 
the wrong Onsager value 1 + h/3 = 17. 4 in the infinite- wavevector limit k — > oo. 

The wavevector-dependent coefficient a T (k) of transverse conductivity has been com- 
puted in two ways, namely, using the relations ^ f£° c T (k, t)dt and ^c 2 k 2 J °° s T (k, t)dt. 
The corresponding results are shown in fig. 3b by open circles and the solid curve, respec- 
tively. These two approaches are mathematically equivalent, but may lead to different 
results in numerical calculations. For instance, the exact infinite-wavelength behaviour 
7/c 2 /47r of function a T (k), where 7 = he 2 lim fc ^ Jo°° S T {^, t)dt, has been reproduced ex- 
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actly by us only in the second approach which, therefore, should be considered as a more 
preferable method of the calculations. As we can see from the figure, the coefficient of 
transverse conductivity a T (k) takes nonzero values for arbitrary nonzero wavevectors, 
contrary to the longitudinal conductivity. 

It is necessary to stress that the transverse dielectric permittivity £ T (k, uj) has a singu- 
larity when both wavevector and frequency go to zero. The asymptotic behaviour near the 
singularity can be obtained on the basis of equation (18) using the equality lim fc ^ £ ' T (k) = 
linn^o e T (k) = e Q . Then one finds that \im k ^o £ T (k, uj) = e Q — vyk 2 /uj. From the last 
expansion it is easy to see that e T (k,uj) is a discontinuous function and its value depends 
on the order of the limiting transitions k, uj — > 0, i.e., lim w ^ h m fc-»o£ T (&;, uj) = e Q but 
lim fc ^ lim^o e T (k, uj) = — ioo. As in the case of longitudinal fluctuations, the singular- 
ity of e T (k,uj) does not lead to singularities in observing quantities and does not break 
the physical requirements imposed on dielectrics. In particular, the macroscopic current 
I^(k,uj) does not appear at uj — > 0, despite the fact that the coefficient a T (k,uj) of the 
polarization conductivity accepts nonzero values at k ^ 0. Indeed, by the definition 
I T (k,uj) = a T (k,uj)E T (k : uj), where the transverse electric field can be defined accord- 
ing to the Maxwell equation via the magnetic field as E T (k,uj) = — j^B(k,uj)xk, so 
that the current vanishes as ~ uj at uj — > 0. If A: tends to zero additionally, we obtain, 
taking into account the asymptotic values 7A; 2 j ' An of a T (k,u), that \im kyLU ^ I T (k, uj) = 
— ^kujB(k)xk J Anc — > 0, where B(k) = lim^o B(k, uj), and the macroscopic current 
goes to zero again without any singularities. 

At small frequencies, the partial contributions of the transverse component of dielec- 
tric functions into observing quantities are small (proportional to uj) with respect to the 
corresponding contributions of the longitudinal part. That is why the functions e T (k) and 
<7 T (k), describing the wavevector dependence of e T (k,uj) and e T (k,uj) in this frequency 
range, have not so important physical meaning as the static longitudinal dielectric per- 
mittivity s L (k). Moreover, since static electric fields are purely longitudinal, the response 
associated with e T (k) and a T (k) cannot be realized phenomenologically in a homogeneous 
isotropic medium, except as the limits 



e T (k) = lim 



E T (k, uj) 



a (k) = lim 



I T (k,uj) 



w-o E Jk,uj) 



We note also that the function e' (k) has no physical meaning for k 7^ at all. This is so 
because lim a ,_ ( .o,fe^o £ ' T (k) / ^(k) = and the wavevector-dependent conductivity a T (k) = 
lim^oi^e^k, uj) — 1) is determined in terms of the imaginary part e'^(k) exclusively. 
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From the other hand, in the infinite-wavelength limit k — > when the conductivity 
vanishes, lim^o <r T (k) = 0, it is not necessary to consider the real part e' (k) as an 
independent quantity because of ]imk^oe' (k) = lim/^o e T (/c) = e Q . 

The normalized time correlation functions $ T (/c, t) = g T (k,t)/g T (k) and T(k,t) = 
s T (k,t)/s T (k), related to dynamical polarization and magnetization fluctuations, are 
plotted in figs. 4 and 5, respectively. The librational oscillations superimposed on the 
exponential, found previously [1, 6] for longitudinal polarization fluctuations, are iden- 
tified for transverse functions $ T (A;, t) as well. But the oscillations damp more quickly 
and their amplitudes are much smaller in this case. The oscillations vanish completely 
at larger wavevectors, namely at k > 20A; min . The magnetization correlation functions 
T T (/c,t) also exhibit oscillatory features, which are more visible at small wavenumbers. 
These functions, however, can accept as positive as well as negative values, contrary to 
the polarization correlations $t(&, t) which remain positive anywhere in time space. It 
is worth remarking that after a sufficiently long period, the polarization functions decay 
purely exponentially in time. This fact has been taken into account by us to calculate 
the contributions of long tails into the dielectric function e T (k, uj) at time integration (14) 
of g T (k,t). In order to demonstrate again that the true choice of microscopic variables 
and AIt f° r describing polarization and magnetization fluctuations in ISM fluids is 
so important, analogous functions, (J 2, T (fc, 0)-J 2 ' T (— k, t)) ^( J 2 T (fc, 0)- J 2 t(— k, 0)) and 
(^#r(fe, 0)-^#r(— k, t)) Q j (^ T (fc, 0)-^#r(— k, 0)) , obtained in the PD approximation, are 
also included in figs. 4 and 5 as dashed curves. The time correlation functions ^(k,t) 
and T(k,t), corresponding to the exact interaction site description (equations (12) and 
(13)), are identical to PD functions in the infinite-wavelength limit k — > 0, but they dif- 
fer between themselves in a characteristic way at greater wavevector values, namely, at 
k > k m \ n . 

The real e' T (k,uj) and imaginary e f ^(k,uj) parts of the transverse dielectric function 
e T (k,u) = e' (k,u) — ie'^{k,u) for the TIP4P water as depending on frequency at fixed 
nonzero wavevectors are shown in figs. 6 and 7 as solid and dashed curves, respectively. We 
mention that lim^o e T (k, uj) = lim^o £ L T (k, uj) = e(uj). For the purpose of comparison 
the corresponding result obtained previously [1] for the longitudinal dielectric permittivity 
e L (k,uj) is also included in fig. 6. The frequency dependence of the transverse dielectric 
function e T (k,uj) at low enough uj can be described by the Debye theory for arbitrary 
wavevectors. This is a result of the fact that the polarization correlation functions g T (k, t) 
damp exponentially at long times. With increasing frequency the collective molecular 
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librations take a prominent role in forming the frequency shape for e T (k, uo), especially at 
small and intermediate wavevector values when k < 20fc min . In this wavenumber region, 
above some uo(k)=l-10 THz corresponding to the position of the first maximum of e'^(k, uo), 
the relaxation behaviour changes into a librational resonance process, characterizing by a 
frequency of order 100 THz. This frequency is associated with the position of the second 
maximum of e'^(k,uo), which practically does not depend on wavevector. At sufficiently 
great frequencies the transverse component tends to the longitudinal dielectric constant, 
so that, for example, at small wavenumbers k < 2/c min the both components become 
indistinguishable from one another at uo > 10 THz (see fig. 6). For larger wavevectors the 
transverse function (fig. 7) differs from the longitudinal one (figs. 6, 7 of [1]) considerably. 
Beginning from frequencies of order uo ~ 1000 THz and wavevectors of order k ~ 100A 1 
the transverse dielectric function tends to the limiting value = 1 corresponding to 
nonpolarizable systems. 

B. Magnetic properties 

The wavevector- and frequency-dependent magnetic susceptibility x(k, uo) = fi(k, uo) — 1 
— —x'(k, uo) — ix"(k, m) has been evaluated for the TIP4P water using fluctuation formula 
(16). Its real x'(k,uo) and imaginary x"(k,uo) parts are plotted in fig. 8 as functions of 
frequency at fixed wavevectors by solid and dashed curves, respectively. We note that 
fluctuation formula (16) is somewhat other by the structure than fluctuation formulas (2) 
and (14) for dielectric functions e L (k,uo) and e T (k,uo). That is why, unlike the dielectric 
susceptibility e L T (k, uo) — 1, the magnetic susceptibility x{k, m) vanishes in the static limit 
at arbitrary wavenumbers, lim^o x(k, uo) = x(k) = 0, and this statement follows directly 
from equation (16). From the physical point of view, such a situation is explained by the 
different nature of electric and magnetic dipoles in ISM fluids. While the electric dipole 
moment di is fixed by the rigid molecular geometry, resulting in = d =const, the 
magnetic moment is caused by the rotational motion of charges sites and it depends 
explicitly on the angular velocity i?j of the molecule. Thus, the vector varies in time 
not only due to the orientational dynamics of the molecule, as d r vector, but also owing to 
the changes of so that |m.j| ^const. As a result, the macroscopic magnetization does 
not appear in the presence of static magnetic fields Hq, because then = even if 
the spatial inhomogeneity of H is taken into account (quasiequilibrium, stationary state 
of the system). 



17 



At nonzero frequencies, when the system is far from equilibrium ((fit) 7^ 0) due to 
the presence of external timely inhomogeneous fields, two mechanisms of appearing the 
macroscopic magnetization are possible. The first mechanism, is connected with the align- 
ment of own magnetic dipole moments along the field Hq, leading to a paramagnetic-like 
behaviour with x'(k, uj) > 0. The second one is caused by the magnetizability of molecules 
owing to the changes of angular velocities in magnetic fields. Then the corresponding ad- 
ditional magnetic moment of the system will be directed oppositely to i3" -vector that 
may lead to a diamagnetic-like behaviour with x'(k, uj) < 0. As we can see from fig. 8, at 
relatively small frequencies, the TIP4P water exhibits paramagnetic features. Beginning 
from frequencies of order 100-300 THz (in dependence on wavevector) and higher the dia- 
magnetic contributions into the magnetization become to dominate and the system under 
consideration behaves like a diamagnetics. It is interesting to remark that the frequencies 
uj p and Ud, corresponding to maximal values of \x'(k, uj)\ in the paramagnetic and diamag- 
netic regions, respectively, almost do not depend on wavevector in a wide wavenumber 
range from k = up to k < 20/c min , where they are equal to uj p ~ 80 THz and ~ 180 
THz, indicating about a weak influence of the translational motions on the processes of 
magnetization. 

Applying the Laplace boundary theorem to fluctuation formula (16) yields the mag- 
netic susceptibility Xoo(k) = h m w^oo x(k, uj) = —hs T (k) < in the infinite-frequency 
regime, where s T (k) = lim^o s T (k, t) > is the static autocorrelation function. The 
function Xoo(^) is plotted in fig. 9 by the solid curve. It takes negative values at arbitrary 
finite wavenumbers and tends to zero as far as k — > oo (note that the magnetic suscepti- 
bility is shown in figs. 8 and 9 for convenience in the negative representation — x)- The 
existence of a nonzero coefficient for the magnetic susceptibility in the infinite-frequency 
regime may be considered as a somewhat unexpected result. For example, the dielectric 
susceptibilities e LT (k, uj) — 1 vanish (as ~ l/^ 2 ) at uj — > oo, because electric dipoles are 
unsensitive to very fast changes of finite electric fields in time owing to the inertness of 
molecules. For the case of magnetic susceptibility the pattern is different for the following 
reason. According to material relations, the macroscopic magnetization of the system 
can be presented in the form j\/l^(k,uj) = ■^-^j^^B^{k^uj). The vector of magnetic 
induction is not independent at k, uj ^ and expressed in terms of the electric field using 
the Maxwell equation as B^(k,uj) = —^kxE(k,ij). So that in the infinite-frequency 
regime the macroscopic magnetization AiT(k,u) vanish as 1/uj at finite values of electric 
fields E(k,u), despite the fact that Xoo(k) ^ 0. A nonvanishing magnetization of the 
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system can be achieved at cu — > oo for infinite values (E(k, uj) ~ a;) of electric fields 
only. But this corresponds rather a hypothetical case which is hard to realize in practice. 
Moreover, it cannot be handled within the linear response theory, where the smallness of 
electromagnetic fields is assumed in advance. 

The value of function Xoo(^) i n the infinite- wavelength limit can be presented in an 
analytical form. Acting in the spirit of derivation of analytical formulas for the static 
current correlations (see Appendix) and taking into account explicit expression (13) for 
the microscopic magnetization density A4t, one obtains the following result 

1 /y- M X,Y,Z /r 1 O -i 1 \ 

Xoo = lim X oo(A:) = -^E^ £ ( J " j~ (K^f + f^A?Af ) , (19) 

a,b a,/3 a a 

where J a are the moments of inertia of the molecule with respect to its three principal 
axes X, Y, Z, the ctth component of vector-position <5" for site a in the molecular principal 
coordinate system is denoted as A" and 1/ J = 1/ Jx + l/ Jy + l/ Jz- It can be shown easily 
that in the PD approximation the function x^P(^) = —h(^ T (k,0)-j£ T (—k,0)) ^2Nd 2 
does not depend on wavevector, i.e., X™(^) — Xoo (the horizontal dashed line in fig. 9). 
The values of x™(^) computed in the MD calculations are shown in the figure as open 
circles. They differ from the exact value Xoo within statistical noise. Using equation 
(19) we obtain for the TIP4P model: Xoo ~ -6.7 • 1CT 10 that is too small m comparison 
with the diamagnetic susceptibility Xh 2 o ~ — 9 ■ 1CT 6 of real water. This is so because 
ISMs do not take into account an electronic subsystem which gives the main contribution 
into the magnetic susceptibility owing to the smallness of the mass m e of electron. This 
contribution can be estimated noticing that the diamagnetic susceptibility of an ideal 
electron gas is defined as Xc ~ —Ne 2 (pl) /Vm c c 2 , where e is the electron charge and (p 2 ) 
denotes the averaged value of square radius- vectors of electrons within the molecule [24]. 
Comparing this result with formula (19) and taking into account that (p 2 ) ~ r 2 , \q a \ ~ |e|, 
max Q [l/J Q ] ~ l/m H r 2 and A a 2 ~ r 2 yields Xoo/Xc ~ m e /m R ~ 10" 4 ~ Xoo/x H2 o' 
where m H is the mass of a hydrogen atom. It is necessary to point out that the atomic 
diamagnetic susceptibility (19) in the infinite-frequency limit, like the case of electronic 
diamagnetism, does not depend on temperature of the system and it is determined by the 
geometry of the molecule. We note that such an atomic diamagnetism can be absent for 
some particular molecular geometries. For example, for the simplest £DS model, when 
M = 2, q 1 = q, q 2 = —q, A 1 = I and A 2 = —I, we find from equation (19) that xf? S = 0. 

We would like also to emphasize that fluctuation formula (16), which explicitly involves 
the microscopic magnetization density, is only the way to calculate the magnetic suscep- 
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tibility x(k, w ) °f ISM fluids in computer experiment not only for nonzero wavevectors 
but also at k — 0. Such a numerical calculation cannot be done within the abbreviated 
description, despite the fact that according to relation (17), the knowledge of the longitudi- 
nal e L (k, uj) and transverse e T (k, uj) components of the generalized dielectric permittivity 
(equations (2) and (6)) allows one, in principle, to determine the frequency-dependent 
magnetic susceptibility x(co>) in the infinite-wavelength limit k — » 0. The reasons for this 
situation are following. First of all we underline that relation (17) is valid for nonzero but 
very small values of wavevector, namely, for k<u/c. Even for relatively great frequencies 
of order 1000 THz this condition corresponds to A;<0.0003A \ At the same time, the 
smallest nonzero value of wavevectors accessible in our simulations is k min ~ 0.3A . So 
that to reach the value 0.0003A 1 it is necessary to increase the size of the simulation 
box to 1000 times! From the other hand, at k = k min we can use relation (17) beginning 
from a frequency of order uj ~ ck m i n ~ 10 6 THz, where e T (k m i n ,u) ~ Xoo(^min)- But it is 
known in advance that the magnetic susceptibility of the system is too small, consisting, 
for instance, approximately —7.5 • 10~ 10 for Xoo(k) at k — k min . Therefore, to reproduce 
this value using the difference of functions e L (/c min ,cj) and e T (k min ,uj) it is required to 
compute them with a relative statistical accuracy of order 10~ 12 at least that constitutes 
an unrealistic problem again. For example, even at extra long simulations with 1 000 000 
time steps, as in our case, the relative statistical uncertainties for the dielectric quantities 
are of order 1%, i.e., 10~ 2 only. 



C. Propagation of electromagnetic waves 

We consider now the question of propagation of free transverse electromagnetic waves, 
E(r,t) ~ Qi(ut-k-r)^ j n foe r jIP4P water. The dispersion relation, connecting frequency 
and wavenumber of the waves in dielectrics, is of the form [18]: 

c 2 -^f^ = 0. (20) 

In this case uj ~ ck and, as it follows from equation (15), the transverse functions e T (k, uj) 
and e T (k,uj) differ from one another by order of x(k,u). As was shown in the preceding 
subsection, the magnetic susceptibility of ISM systems is sufficiently small at arbitrary 
wavevectors and frequencies, i.e., x(/c, cj) <C 1. Thus we can putt e T (k,uj) ~ e T (k,u) 
in relation (20) without loss of precision. Further, the characteristic frequency scale of 
varying the dielectric constant is of order 1000 THz. This corresponds to an interval of 
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varying wavenumbers in the waves of order 0.0003A . At the same time, the character- 
istic wavevector scale of varying e T (k,uj) is of order k min ~ 0.3A \ so that the spatial 
dispersion of the dielectric constant in equation waves (20) can be neglected completely, 
i.e., e T (k,u) ~ e T (k,ou) ~ e(u). The dielectric permittivity e(u) is an imaginary function 
of frequency. Therefore, a solution to the equation (20) with respect to wavenumbers at 
a given real frequency (it corresponds to frequency of the source which creates electro- 
magnetic waves) it is necessary to find in the form k = k' — ifc", where the magnitude k" 
of the imaginary part will describe the damping of waves with the wavelength A = 2ir/k' . 
We shall consider a case when plane waves are damped in direction of their propagation, 
i.e., when fc'-fc" = k'k". 

Numerical results for the dispersion u(k') and damping k"(u), coefficients in the in- 
frared region of spectrum are shown in figs. 10a and 10b, respectively. We note that 
calculation of the damping coefficient at great frequencies is very sensitive to the preci- 
sion of evaluation of the dielectric constant. A more accurate estimation can be achieved 
when the dielectric constant is determined in terms of current fluctuations via fluctuation 
formula (5) (the solid curve in fig. 10b), instead of the usual formula (2) (the open circles). 
As far as k" ^ in the whole region of frequencies < uj < oo, purely monochromatic 
waves cannot propagate in dielectrics. Nevertheless, choosing a criterion k" <C k', we 
can talk about quasimonochromatic waves with a slight absorption. According to our 
calculation, only two opposite regions of very small (cj^O.01 THz; radiowaves) and very 
great (u ^ 500 THz; far infrared, visible light) frequencies can satisfy this criterion. In 
the radiowaves region (see insets of the figures), the phase velocity uj(k')/k' = c/ \feo is 
defined by the static dielectric constant ^fe^ pa 7 (for real water y/e^ pa 9). For example, 
at uj w 0.01 THz we obtain k'c « 0.1 THz and k"c w 0.002 THz. Thus, the radiowaves 
with the wavelength A = 2ir/k' pa 2cm are dumped in a characteristic way during the 
interval 2n/k" pa lm. In the far infrared region the phase velocity of the nonpolarizable 
TIP4P water tends to its infinite frequency value c/y 7 ^, where = lim s(uu) = 1. 
Owing to the electronic polarizability of molecules, which is not taken into account in the 
TIP4P model, the dielectric constant of real water differs from unity even in the visible 
light spectrum uj ~ 3000 THz, where \Js(uj) ~ 1.33. It is interesting to remark about the 
existence of a large region of frequencies (uj ~ 10 — 50 THz), where the group velocity 
duj{k')/dk' = c/\fe* is practically constant. In this region, near the librational minimum 
of the imaginary part e"(u), the real part e'{uS) of the dielectric constant approximately 
does not depend on frequency and takes values of order e* pa 2.3 (see fig. 6a of Ref. [1] 
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for s{uj) and fig. 6a of the present paper for e T (k m i n ,uj) ~ £(<^))- However, the waves are 
damped significantly (k" < k') in this frequency range. Within the absorption maximum, 
which occurs at ~ 150 THz, we have k" ~ k', so that electromagnetic waves decay here 
on an interval of order of their wavelengths. 

All our previous calculations of the dielectric functions e L T (k, uj) dealt with real values 
of wavevector and frequency. Such functions at a given single set of k and uj describe the 
response of the system on electromagnetic fields in the form of monochromatic plane 
waves ~ e i(wt-fc-r)_ imaginary values of frequency and wavevector correspond to cases 
when amplitudes of fields either increase or decrease in time and space. Since arbitrary 
inhomogeneous fields can be cast as a set of the monochromatic waves, using the time and 
spatial Fourier transform, the dielectric function at imaginary values of k and uj may be 
expressed in terms of its values at real wavevectors and frequencies. As a demonstration, 
we consider the case of purely imaginary frequencies, uj = icj*, in the infinite wavelength 
limit (k — 0). Then, as it follows from fluctuation formula (2), the dielectric constant can 
be calculated directly 



S(t)di , (21) 




using the time correlation function $(£) = {J2ij di(0)-dj(t)} / di(0)-dj(0)) of the 
total dipole moment. From the other hand, presenting the field ~ e~ w ** in the form of 
the corresponding Fourier integral over real frequencies, it can be shown, in particular, 
that at uj* < the dielectric constant e(u*) is expressed via the imaginary part of e(uj) 
(defined at real frequencies) as follows [17]: 



00 



s(cu*) - 1 = - / X f iX \ dx . (22) 
nJ x 2 + uj* 2 





The function e(uj*) is plotted in fig. 11a. It describes the response of the system on 
external electric fields which increase in time exponentially. As was expected, this function 
at uj* — > —00 tends to unity owing the inertness of polar molecules. 

The case uj* > will correspond to decaying in time of electric fields. This decaying 
can be caused either by decreasing of external fields or by switching off of sources which 
support these fields. As an example, consider the spatial and time decay of electromagnetic 
fields in the system after the passage of transverse electromagnetic waves. In this situation, 
frequency and wavevector are purely imaginary quantities, i.e., uj = iuj* with uj* > 
and k = —ik", that is necessary to take into account finding a solution to the equation 
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(20). For the same reasons as in the case of propagation of electromagnetic waves, we put 
e T (k, uj) e{uj) in equation (20) and obtain uj* = ck" / \Je(u*), where e(uj*) accepts purely 
real values at purely imaginary frequencies. We note that an independent parameter in 
this equality is k", so that the function oj*(k") will describe the time decay ~ e~ w ** of 
the transverse electric field with a spatial inhomogeneity of k" which has been created 
by the passed electromagnetic wave (see fig. 10b). The function uj*(k") is shown in 
fig. lib. As can be seen, weak spatial inhomogeneities correspond to long life times 
~ l/u*(k"), whereas strong inhomogeneous fields damp in time faster. It is worth to 
remark also that in the asymptotic limit t — > oo, the time correlation function $(t) decays 
in time exponentially as ~ e~'/ rro1 , where r rc i = 6.7ps is the relaxation time. Therefore, 
the dielectric constant e(u>) at imaginary frequencies can be defined only in the region 
^suj < l/r re i = 0.149 THz, because otherwise the integral in (21) is divergent. This merely 
means that transverse electromagnetic excitations cannot damp in time faster than with 
the characteristic interval r re i. The limiting region of imaginary frequencies is shown in 
fig. lib by the horizontal dashed line. 

Finally, as far as real e' and imaginary e" parts of the transverse wavevector- and 
frequency-dependent dielectric constant are defined according to the fluctuation formula 
(6) via the same time correlation function c T (k,t), they are not independent and must 
be connected between themselves. The desired relations can be obtained using analytical 
properties of the functions c T (k,t), e' T (k,uj) and e'^{k,uj) — Aita^k) / uj. The result is 

_ 1 = I l x ^4 &x , 4M = *± /^^HW . 

1 TV J X — UJ 1 7T J X — Ul UJ 



(23) 

The relations (23) are similar to the well-known Kramers- Kronig expressions for the di- 
electric constant of conductors [17]. They can be applied to the transverse frequency- 
dependent dielectric constant of ISMs at arbitrary values of wavevector. 



IV. CONCLUSION 



We have established that the calculation of the dielectric quantities for interaction 
site models of polar fluids can be performed in computer experiment by two alternative 
ways, namely, using either charge or current fluctuations. The first way is more efficient 
to evaluate the longitudinal component e L (k,uj) of the dielectric permittivity, whereas 
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the second method is suitable to obtain its transverse part e T (k,u). Separating the 
total current of moving charges into electrical and magnetic components, the transverse 
generalized dielectric permittivity e T (k, ui) = e T (k, u) + ~l~ l+^i) nas been expressed in 
terms of the transverse dielectric e T (k, u>) and magnetic x(k, oS) functions. These functions 
have been evaluated for the TIP4P model of water by molecular dynamics simulations 
in the whole region of wavevectors and frequencies for the first time on the basis of the 
proposed fluctuation formulas. 

It has been shown that the transverse e T (k,u) and longitudinal e L (k,u) = e L (k,u) 
components of the dielectric tensor e(k, uj) differ between themselves in a characteristic 
way and coincide with one another in opposite limits of either very small wavenumbers or 
very large wavevector and frequency values. Moreover, contrary to the case of longitudinal 
fluctuations, the generalized dielectric permittivity e T (k,u) exhibits a specific behaviour 
when both wavenumber and frequency tend to zero. We has identified also that at great 
frequencies the TIP4P water can be considered as a weak diamagnetics. However, values 
of the magnetic susceptibility x(k,u) are much smaller in amplitude than those for real 
water, because the electronic magnetization is not included in this model. It is worth to 
stress that in the present investigation the polarization and magnetization microscopic 
densities have been constructed taking into account explicitly the atomic structure of 
ISMs. As a result, it has been demonstrated, in particular, that at nonzero wavevectors the 
genuine transverse dielectric function e T (k, u) has nothing to do with that obtained in the 
point dipole approximation. The last function behaves like the dielectric permittivity of a 
Stockmayer system [12-14] and can be used at small wavenumbers only as an estimation 
of the frequency-dependent dielectric constant in the infinite-wavelength limit. 

The longitudinal and transverse components of the wavevector- and frequency-depen- 
dent dielectric permittivity describe all electromagnetic phenomena in the system. The 
knowledge of these quantities may present an interest in both theory and pure experiment. 
The performed calculations can be extended, in principle, to more realistic models of polar 
fluids. This will be the subject of a separate consideration. 

The author would like to acknowledge financial support by the President of Ukraine. 
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Appendix 



We shall show that the longitudinal and transverse components c L T (k) of the static wavevec- 
tor-depended Kirkwood factor of second order c(k) = c L (k) + 2c T (k) are presented analytically. 
According to definitions (3) and (7), the functions c{k) and c L (k) read 

1 lN,M N,M \ 

\ i,a j,b I g 

(Al) 

1 lN,M N,M \ 

\ i,a j,b I q 

where the site velocities can be split as V" = Vi + f2iXSf. We mention that Vi and fii are 
the translational and angular velocities of the ith molecule, respectively, 81 = rf — Vi and r« 
denotes the position of the molecular centre of mass. 

Equilibrium distribution functions are factored into the coordinate and velocity parts. In 
its turn, translational and angular velocities are distributed independently of one another for 
each molecule. As a result, nonzero contributions to equilibrium averaging give only terms with 
coincident molecular indexes {i = j) of summation (Al). It is more convenient to consider each 
molecule in its own principal coordinate system in which 8^ = A a . Then expressions (Al) 
transform into 

-, N M 

i a,b 

(A2) 

-, N M 

c L« = ]v^EE9A<(( fe ^0 2 + (fc-[^xA a ])(fc.[^xA 6 ]))e- ifc -^) o , 

i a,b 

where p ab = A a — Aj. It is necessary to underline that for a given molecular geometry, A a 
(a = 1,...,M) constitute the set of constant vectors characterizing the positions of charged 
sites in the principal coordinate system attached to the molecule. Further, we use the follow- 
ing equalities 3(k-[f2x A a ])(fc-[17x A b ]) = (3kk - 1) : fix A a 17x A b + [17 X A a ] • [fi X A b ] and 
[fix A a ]-[i?x A b ] = f2 2 (A a -A b ) - (fi-A a )(fi-A b ) and take into account the relations 

( e " fc ' P )fc = jo{kp) ' ( (3 ^ " 1)e " fc ' P )fc = -^ k P^PP ~ 1) , (A3) 

where j Q (z) = sin(z)/z, j 2 (z) = Sj^z)/ z — j Q (z) are the spherical Bessel functions of order 
zero and two, respectively, p = p/p and averaging in (A3) is performed over orientations of 
fc-vector with respect to the molecule. Owing identity of molecules and isotropy of the system 
we have, in particular, (V?) = (V 2 ), (f2?) = {Q 2 ) and ((k)-V) 2 ) = (V 2 )/3. Using the definition 
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of temperature and the equipartition theorem yields m(F 2 )o/3 = J a (^a)o = k B T, where m is 
the mass of the molecule and Q a are the principal components (a = X, Y, Z) of angular velocity. 
Finally, in view of the statistical independence of angular velocities directed along different 
principal axes of inertia, after cumbersome but not complicate operations one obtains 



k B T M 

a,b 



— V (- - — 



(A4) 



C L {k) = ^c(fc) - ^ J2 1a%j 2 (kpab) Kb ( 3 PabPab ~ <W) > 
aj^b a,f3 

c(k) - c L (k) 1 k B T^ f 

C T {k) = ~ = -C(fc) + 2^ qaQ b J 2 {kpab) 2^ Kb [ZpabPab ~ <W 

aj^b a,f3 



(A5) 



where = -AfA£/J 7 and h%? = Af Af /J 7 + A^Aj/J^ (all the three variables a, (3 and 7 
take different values here) denote the constant quantities which are defined by the molecular 
geometry and p a b is the distance between sites a and b within the same molecule. 

Thus, the static Kirkwood factor of second order is uniquely determined by the temperature 
and geometry of molecules. It is the same for infinite and finite systems, i.e, Cl,t(&) = c L T (^) 
(infiniteness of the system has not been used at the derivation of equations (A4) and (A5)). As 
can been shown earlier, this statement is completely in line with transformations (9). 
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Figure captions 



Fig. 1. The wavevect or- dependent longitudinal Kirkwood factor g L (k) (a) and the 
longitudinal and transverse components c LT (k) (b) of the current autocorrelation func- 
tion for the TIP4P water. The dashed and solid curves in (a) correspond to the finite 
and infinite systems in the ISRF geometry. The results of the Ewald geometry for the 
longitudinal and transverse components are presented by the full squares and circles, re- 
spectively, whereas the corresponding data in (b) obtained within the ISRF geometry 
are displayed as open squares and circles. The analytical evaluation of c LT (k) (equation 
(A5)) is plotted in (b) by the solid curves. The calculations of g L (k), performed through 
the dynamical current correlations (see equation (8)), are shown by the direct and rotated 
crosses for the Ewald and ISRF geometries, respectively. 

Fig. 2. The normalized dynamical total current correlation functions of the TIP4P 
water. The result of the Ewald geometry for the longitudinal and transverse components 
is presented as open and full circles, respectively. The data, obtained within the ISRF 
geometry for the infinite system, are plotted by the corresponding solid curves. The 
longitudinal component of the finite system within the ISRF geometry is shown as the 
dashed curve. 

Fig. 3. The wavevect or- dependent transverse dielectric function e T (k) (a) and gener- 
alized conductivity a T (k) (b) in the low frequency limit (the solid curves). The dielectric 
function corresponding to the PD approximation and the conductivity, calculated in terms 
of total current (instead of magnetization) correlations are shown in subsets (a) and (b) 
as the dashed curve and open circles, respectively. 

Fig. 4. The normalized, time autocorrelation functions of the transverse polarization 
fluctuations in the TIP4P water. The results of the PD approximation are shown as 
dashed curves. Note that in the infinite- wavelength limit k — > (see subset (a)), the PD 
functions are identical to genuine ones. 

Fig. 5. The normalized, time autocorrelation functions of the transverse magnetiza- 
tion fluctuations in the TIP4P water. Notations as for fig. 4. 
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Fig. 6. The transverse component of the wavevector- and frequency-dependent dielec- 
tric function e T (k,cu) of the TIP4P water at small wavenumbers. The real and imaginary 
parts are plotted by the bold solid and dashed curves, respectively. For comparison the 
longitudinal component is shown by the corresponding thin solid and dashed curves. 

Fig. 7. The transverse component of the wavevector- and frequency-dependent di- 
electric function of the TIP4P water at great wavenumbers. The real and imaginary parts 
are plotted by the solid and dashed curves, respectively. 

Fig. 8. The wavevector- and frequency- dependent magnetic susceptibility of the 
TIP4P water. The real and imaginary parts are plotted by the solid and dashed curves, 
respectively. Note that the results are shown in a negative representation. 

Fig. 9. The wavevector- dependent magnetic susceptibility in the infinite-frequency 
regime (full circles connected by the solid curve). Note that the result in the PD ap- 
proximation is independent on wavenumbers (horizontal dashed curve). The PD function 
calculated in the MD simulations is shown as open circles. 

Fig. 10. The dispersion (a) and spatial damping (b) of transverse electromagnetic 
waves in the TIP4P water. 

Fig. 11. (a) The dielectric constant of the TIP4P water in the infinite- wavelength 
limit as a function of purely imaginary frequencies, (b) The time decay of electromagnetic 
fields in the system after the passage of transverse electromagnetic waves. 
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